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The kinetic theory of a self-gravitating system is considered in the Bhatnager- 
Gross-Krook approximation to the kinetic equation. This approach offers a unique 
and tractable setup for studying the central, collision-dominated region of the 
system, as well as its almost coUisionless outer part. Stationary non-equilibrium 
states of the system are considered and several different self-similar solutions are 
identified. 



1 Introduction 

There is a long traditipn of applications the kinetic theory to the gas of self- 
gravitating particles tl-Q. The first regular approach of that kind was developed 
by Chandrasekhar El, and appeared to be very useful in astrophysical applications. 
Indeed, the kinetic theory is certainly applicable for globular star clusters, which 
are rather well-isolated systems on the relevant time-scales, and contain sufficiently 
many (typically of order N ~ 10^) stars a. Therefore, the kinetic approach, which 
typically studies a closed, macroscopic system, has a definite range of validity here. 
It allows to identify correctly the coUisional time-scalexl-Q, and predict the moderate 
and late stages of the evolution for globular clusters 

The general method for constructing the coUisional Fakker-Planck equation for 
statistical systenis with 1/r interaction was proposed in El (see also the text-book 
descriptions in Elaa^Dfl). This equation adopts a more general strategy for Boltz- 
mann's kinetic equation to the considered long-range interacting system. After 
suitably dealing with the peculiarities of the long-range interaction, one ends up 
with the coUisional Fokker-Planck equation. The main difficulty with this equation 
is that it is technically rather involved and resistant to reliable analytical treat- 
ments. Thus, one has to change the level of description and move to hydrodynam- 
ical models of the stellar dynamics (see e.g., Q'Q and refs. therein). Instead of the 
one-particle distribution function, which is the basic object of the kinetic theory, 
the hydrodynamical approach operates with a few coarse-grained hydrodynamic 
variables, but by its meaning it is applicable only to the dense central regions of 
the self-gravitating system, where collisions are certainly dominating. 

In this contribution we present another approach to the kinetic theory of self- 
gravitating systems. This is an approach of kinetic models, where instead of going 
immediately to the hydrodynamical level of description, one is properly modelling 
the coUisional part of the kinetic equation, trying to preserve its main qualitative 
characteristics, but to make the system analytically tractable. One of the most 
Pj^LCps^afnl kinetic models was proposed by Bhatnager, Gross and Krook (BGK) 
f"^^' and since that time appeared to be very useful in the kinetic theory. 
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The reason of its success is that the model offers a niinimal way to incorporate 
the necessary conservation laws (those of particle number, momentum and en- 
ergy) with a simple relaxation mechanism. For a self-gravitating statistical system 
this approach provides natural methods to consider simultaneously the collision- 
dominated (overdamped) regions of the system, and its almost collisionlcss outer 
domains (underdamping) . 

In the present contribution we will restrict ourselves to the presentation of the 
BGK-method as applied to self-gravitating systems (section II), and to straight- 
forward analytical investigations which will recover within a single setup several 
known and also new results for the steady regime of behavior (section III). More 
detailed investigation of the situation, as well as the presentation of more elaborated 
BGK-type models, are planned for future. 

2 Bhatnager-Gross-Krook kinetic equation 

The statistical dynamics of N point particles with unit mass interacting through 
Newton's inverse-square law, will be considered in the classical kinetic approach. 
This basically involves two important assumptions: (i) The field acting on a given 
particle (test particle) can be represented as a mean, self-consistent field plus a 
contribution from two-particle collisions, (ii) As usual, a collision between a pair 
of particles is taken to be independent of the others, and local (namely, it leads 
to sudden changes in the corresponding momenta, whereas the coordinates can be 
considered as fixed). 

The mean-field assumption is reasonable for a system having long-range forces, 
since they suppress fluctuations. On the other hand, the second assumption con- 
cerning the simple cumulative effect of independent collisions was subjected to 
a certain criticism (see e.g., 0), because the long-range unshielded gravitational 
coupling involve many simultaneously interacting particles. Nevertheless, an im- 
pressive amount of observational and numerical material has been obtainedL which 
advocates the use of the present approach al least on certain time-scales uij. 

According to our assumptions the contribution from the direct interaction be- 
tween the particles is added on the right-hand side of the Liouville equation for the 
one-particle distribution function /(r, v, t) as the corresponding collision integral 



The mass of each particle is taken be unity, m — \, and the potential (/)(r) can be 
viewed as an external despite its self-consistent closure with the Poison equation: 



As usual in kinetic theory, the distribution function /(r, v, t) is normalized on the 
total number of particles, namely 




(1) 





(3) 
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is the density of particles, and 

J drn{r) = N (4) 

is the total number of particles. 

Instead of deriving the collision integral of the kinetic equation from the under- 
lying microscopic dynamics, Bhatnager, Gross and Krook proposed the following 
model for the coUisional integral □: 

CBGKif] = -iy{r)[f{f,v,t) - /o(r,i/,t)], (5) 

where iy{r) > 0, and fo is the best local Maxwellian 

where n{r, t), T[r, t), u{f, t) are defined by the following conditions: 

n{f,t) - / &v f{r,v,t), (7) 



{r,t)u{r,t) = / dvv f{r,v,t), (8) 



3n(f, i)T(f, i) = j dv{v~- uif, t)Y fir, v, t). (9) 

Due to this conditions the collision integral conserves probability, momentum 
and energy, and in the spatially homogeneous case T = T{r, t) it enforces conver- 
gence of the distribution function /(r, v, t) to the corresponding Gibbs distribution. 
One should notice that BGK collision integral (||) is still non-linear due to these 
self-consistency relations. However, this non-linearity is based on the first three mo- 
ments of / only, which means that BGK scheme will be friendly to moments equa- 
tions and related procedures. Another attractive property of BGK model is that 
though being proposed ad hoc, it was later demied from the Boltzmann equation in 
both rarefied (low-density), and dense limits Ei£3. In particular, relations with more 
formal approximatkm methods (e.g., those proposed by Chapman-Enskog or Grad) 
were established Ii2t2l. Thus, there is a general expectation that this approximation 
will reasonably mimic the main properties of the original Boltzmann equation. 

The general conditions of invariance and convergence do not specify v{r), which 
has a meaning of the inverse characteristic relaxation time. It should not depend 
on V, since else the conservation laws cannot be satisfied. Obviously, the concrete 
choice of v{r, t) should be determined demanding that Eq. (||) is as close as possible 
to the original Boltzmann equation. A standard assumption is that ly depends on 
the coordinate r and time t through the r-dependence of the first three moments, 
: v{n,u,T), and a general procedure was proposed to derive its concrete form 
We will not go into details of that procedure, but we will show that already 
rather simple qualitative considerations determine its rough form, and then we will 
directly take for it the expression for the inverse diaracteristic time, which was 
obtained from the full Boltzmann equation in Refs. ui3. To determine v one notices 



extended'version: submitted to World Scientific on February 1, 2008 3 



that in the rarefied case it is proportional to n, and then due to obvious dimensional 
reasons one has 



v{r, t) cx 



n(r, t) 



(10) 



r(f,t)3/2' 

In this context one could wonder why only T and not <j) should enter v. Since u 
belongs to the collision integral, according to the common initial assumption the 
influence of the field is not taken into account there. Indeed, in the full Boltzmann 
equation binary collisions are considered in such a way that only microscopic two- 
particle interaction enters there, and not the mean-field. This is inevitable to find 
tractable schemes, and means that collisions occur on such time and length-scales, 
where the effect of the mean- field is not relevant. 

More precisely, we notice the bmajcv. collision relaxation time obtained through 
the complete Boltzmann equation Gl'Ela'Q 

0.065(v2)3/2 



11 



nraG^ ln(0.4iV) 

with the right hand side calculated at typical radius R, and where we momentarily 
recovered the dimensional units. Thus, we just take 1/v equal to the local relaxation 



time, so 



v{r) 



1 



nmG^ ln(0.47V) _ 15.4 ln(0.47V) G'^nt'/'^nir) 



7 



tr{r) 0.065(t;2)3/2 
15.4 ln(0.47V)G^TO^/2 



r(r)3/2 



nir) 



(13) 



is the inverse relaxation time at position r a. To estimate whether collisions are 
relevant compared with the motion in the mean-field, one has to compare i^ei with 
the dynamical time-scale td- For a cluster with mass M = Nm and typical radius 
R the dynamical time is 

i?3/2 ^ ^ 

(14) 



_ R _ 



The situation tr{r) <C td means that collisions dominate (overdamping), and 
tr{r) ^ td indicates underdamping. This qualitative criterion will be applied to 
determine the transition from the outer weakly-damped part of the system (halo) 
to the inner strongly-damped part (core). 



2.1 H-theorem and convergence to equilibrium 

Let us now briefly show that if the conditions for equilibrium are satisfied, then the 
BGK-collisional term (||) indeed induces convergence toward the Gibbs distribution. 
To this end we will consider the time-behavior of the Boltzmann-Gibbs entropy: 



S{t) 



dr dv / (r, t/, t ) In / (r , w , t) , 



(15) 



S 



dr dw z^(r)[/ — /o] In/ + / dr dv 



d(t)_dl 
dr dv 



In/ 



(16) 
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Using the fact that J dv [f — /o]ln/o = 0, which is conservation of probabihty, 
momentum and energy, one can rewrite the first term in r.h.s. of Eq. ([l^ ) as 



/hi^ + /oln^ 
Jo J 



>0, 



since 



dv f\n^ > 
/o - 



(17) 



(18) 



for any distributions /, /o- The second term term in r.h.s. of Eq. (^^ can be 
presented as an integral over the corresponding surface in the coordinate space 



dS / dv vf In / 



(19) 

If there are no global currents this integral is zero, and we obtain the H-theorem: 

S>0 (20) 

Monotonically increasing at finite times S will saturate in the infinite-time limit: 
f = fo- Substituting this expression in l.h.s of the kinetic equation we will get the 
Gibbs distribution 



T = const. 



n{r) = 
^ ^ Z 



(21) 



If there are global currents in the system, then the actual distribution is not known 
a priori, and the full kinetic equation has to be solved. 



2.2 Kinetic equation in spherical coordinates 

Since we will consider non-equilibrium distribution functions, which are in a sense 
close to be isotropic, it is reasonable to employ spherical coordinates: 



X = r smtl cos ip, 



y — r sma smip, 



z = r cos f. 



(22) 

In the local Cartesian base {eV, ee, e^p} (e.g., eV is a unit vector, which points out 
from r to r -|- dr) velocity will have components 



where 



V = VrCr + vgee + v^e^, 



Vx = Vr Sm U cos if + Vg cos U cos if — V^ sm if 



Vy = Vr sm t) sm Lp + vg cos sm Lp + v^p cos ip, 
Vz — Vr cos 9 — Vg sin 0, 

In the spherical coordinates the kinetic equation will read 



(23) 

(24) 
(25) 
(26) 



dl 
dt 



-Vr 



df , Vg df 



V, 



df 



dr 
cote 



r do 

dvg 



r sin 9 dip 



Vn 



V, 



df_ 

dVr 



V^pVg 



dl_ 

dv,„ 



Vr 

r 



Vg 



dl 

dvg 



dl 

' dv,„ 



= CBGK[f] (27) 
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We will consider non-equilibrium states with radial symmetry. We will assume that 
(j}{r,9,ip) = </'('"), and the distribution function / is non-equilibrium exclusively 
through its radial properties: 



/(r, e, if, Vr, ve,v^) = f{r, Vr) 



1 



27rr(r) 



cxp 



2T{r) 



(28) 



namely there is no dependence on the angular coordinates, and the angular veloc- 
ities have relaxed to equilibrium. In fact, Eq. ( p8| ) is the minimal Ansatz, which 
allows to study non-equilibrium states within the simplest spherical geometry. 
For the collision term one has 



CsGKif] = -Hr) [f{r,Vr) - fo{r,Vr)] X 



27rT{r] 



exp 



' 2T(r) 



(29) 



foir,v) 



n(r) 



[27rr(r)]i/; 



■ exp 



{v — u)' 
" 2r(r) 



(30) 



where the mean radial velocity u describes an average motion (current of particles) 
along the radial direction. 

Let us now substitute Eq. (|2|) to Eq. {^t\), and integrate by vg, v^. The result 
will read 



df ^ df 2T(r) 

Ot OVr r 



1 (r) OVr 



^v{r)[f(r,vr)~ h{r,vr)] (31) 



where /' = df{r,Vr)/dr, (p' = d(f>/dr. 



2.3 Moments 

Eq. ( ^ ) is still complicated integro-differential equation. To subtract the relevant 
information from it one considers the subtracted and non-subtracted moments 



Mi{r) = / dVr {Vr ~ uf f{r,Vr), 



M 



(r) = J dvr vl.f{r,Vr), 



(32) 



(33) 



These two types of moments correspond to different experimental situations, since 
Mi{r) are measured in the comoving frame, whereas Mi at the rest. We will use 
them in parallel choosing the one or the another as appears most convenient. 
The corresponding equations read: 



1 



Ml + liiMi^i + — [r2(Af,+i + uMi)]' + A/,_i + lu'[Mi + uM,_i] 



2T 



^ -v{r)[Mi - uj{l) nT"^ {I - l)% 
where uj[l) is zero for I odd, and equal to 1 for / even. 



(34) 
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1 / 2T\ 
Ml + -{r^Mi+iY + M 0' Mi-i 

I 

= -i,{r)[Mi - nY^u^-PTP/^ u;{p){p-l)U], (35) 

p=0 

Recall that, by construction, the low moments 

Mo^n{r), Mi = 0, M2 = n(r)r(r), (36) 

Mo = n{r), Mi=n{r)u{r), M2 = n{r)[T + u^], (37) 

are the same for fo{r,Vr) and f{r,Vr). Since this fact expresses conservation of 
probability, momentum and energy, the equations (|3j, ^ ) with I = 0, 1,2 will be 
universal, namely they do not depend on the concrete form of the collision integral. 
Notice also that the third moments are connected as: 

M3 = A/3 + nu^ + 3unT. (38) 

Here A^3 is energy current, which consists of heat current M3, transfer of kinetic 
energy nuu'^, and work done by pressure 3>unT = iuP. 
Let us write down the first members of Eq. (jsj) 

/ = 0: r2n(r) + (r^nw)' = 0, (39) 
1 = 1: nil + (nT)' + n(j)' + iiuu = 0, (40) 

1 = 2: nt + ^(r^MsY + 2u'nT + unT' = 0, (41) 
1 = 3: M3 + inTii + \[r^{Mi + ulVh)]' + iu^Ah + nuT] 



+3nT{4>' U) = _^(^)M3, (42) 

r 

l = A: M4 + 4M3U + [r^ (Ms + uA'U)]' + 4u'[M4 + uM^] 

+4M3(</>' ^) = -u{r)iA'h - 3nT^) (43) 

r 

Eq. (|39| ) expresses conservation of the mass inside the sphere r. When u = 
Eq. (40) expresses hydrostatic equilibrium; when it 7^ it contains flow terms, that 
should be easy to understand. The equation for I = 2 expresses the conservation of 
energy; for obtaining it we have also inserted Eq. (|39|). 

Somewhat different, but of course, totally equivalent expressions will be obtained 
for the Al-moments. 

Eqs. ( ^9|^3| ) are inherently attached to the Poisson equation: 

^{r^(j)'y = 4:TTGn{r). (44) 
The density n{r) should satisfy to the condition of integrability, 

/•OC 

47r / drn(r)r^ = N. (45) 
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Notice that equations with I > 2 can be written in a form, which does not 
contain (j)' exphcitly. In that respect their combination can be viewed as equations 
of state. For example, using Eq. po|), the for result for / = 3 reads 



1 



Ah + :^[r {Mi + uAh - SnT^)]' + Su'JVh + inTT' = -v{r)M:i. (46) 
3 Steady state. 

As was already indicated in the introduction, we will consider the stationary 
(steady) situation. This case, where all quantities are time- independent, is worth 
to study for its own sake, since there are clear observational evidences that typical 
globular clusters do have such a phase in their evolution I3'l3. On the other hand, it is 
a necessary step towards understanding of more general, time-dependent situation. 

In the situation, where the known conditions for equilibrium are satisfied (e.g., 
absence of energy and/or probability currents), the self-gravitating system under 
study has the Gibbs distribution as its only steady state. This point was already 



illustrated by us when considering the H-theorem in section 2.1. In certain range 
of temperatures this distribution is physical (e.g., it is at least metastable)_but it 
loses its stability for temperatures lower than a certain critical tempeijatureQ. This 
is the notorious phenomenon of gravo-thermal instability (collapse) EIqiO. However, 
there are situations, where one expects that the very reasons for the existence of 
the equilibrium can be invalid. The most typical situation of that kind appears 
with binajY star formation and tightening processes in the central region of a star 
cluster Bq'EI. These processes are accompanied with a release of energy, so that 
the rest of the system (namely its outer with respect to a relatively small part, 
where the binary-formation process occurs) can be considered as being subjected 
to sources of energy put in the central part. The existence of energy currents from 
the central part leads to a non-equilibrium steady state. Additionally, they can 
stabilize the behavior of the self-gravitating system preventing its collapse, since 
the mechanism of the gravothjecmal instability is connected with spontaneous energy 
currents towards the center Ij^Bq. 

A somewhat more exotic example of a non-equilibrium steady state cap 
provided by the existence of a black hole in the central part of a star cluster 
Then one has a stationary current of consumed stars, which is directed towards the 
center, as well as a related current of energy. 

Our setup with the moments equations is especially suitable for describing the 
above non-equilibrium states, since there are explicit expressions for the energy 
current, M^^ and the current of particles u. 



3.1 Ideal hydrodynamics 

This scheme is determined by a condition M3 — (no heat current). The nonequi- 
librium character is displayed only through the stationary current of particles u. 
Eqs. (p9l-P) read 



(r^n u)' = 0, (nT)' + n0' + nuu' = 0, 2u'nT + unT' = 0. (47) 
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Recall that the usual gibbsian equilibrium is a particular case of this situation, 
which is realized for m = 0, and an additional assumption that T — const. In this 
case Eqs. are trivial, and Eq. (|4^) leads to the Gibbs distribution. 

Let us first assume that u does not depend of r. Then Eq. (|4^) predicts that 
T = const. With Eq. @ one gets 

c 1 , 

n= 7T, c = const, (48) 

u 

Notice that similar to the isothermal case this distribution is non-integrablc at 
infinity, and therefore predicts infinite mass. Being combined with the Poisson 
equation (|4^), Eq. (|4^) offers an exact solution: 

= 2rinr + (j!)o. (49) 



To study the case u ^ const one first notice that Eq. (47) can be exactly 
integrated 

T^%, ci= const (50) 

Eqs. (H, |ol) combined with (|^]4|) lead to 

ciur^{u~^r-^Y(t>' + u'u = 0, u{r^(j)'Y ^ iirGc. (51) 

To study self-similar solutions, one makes an Ansatz u = uor", and gets from 
Eq. 

-ciw-2(3a + 2)r-2"-i +0' + Q,u2^2a-i (52) 

Consulting with the Poisson equation one observes that for small distances the last 
term in Eq. ( [5^ ) can be neglected, and as the result one has the following asymptotic 
solution at small distances a = 1/2 and 



7ci _2 , SttGc _3/2 Ci _i C _5/2 

2ui uo Un uo 



It is seen that the density is singular at the origin, but still integrable. The first 
term in 0' indicates the presence of the central mass 7ci/{2uqG). This implies 
the following choice: uq < 0, c < 0, and then the obtained solution describes a 
stationary consumption regime of stars by the central mass. The energy current 

3 



Ads, = nu^ + 3nuT is also directed towards the center. 



3.2 Overdamped situation 

Starting since this moment, we will into the consideration an energy current, de- 
scribed by M3. We make a self-consistent assumption that there is a range of 
distances close to the inner part (i.e., the core), where collisions dominate. In this 
overdamped regime 7 7 = (dynamic time) /(kinetic time) is large, and the large iy{r) 
in the r.h.s. of Eqs. ( [42| , [43| ) imposes the Gaussian behavior for the corresponding 
moments: 

M3 = 0(1/77), M4-3nT2 = 0(1/77) (54) 
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Using the last result as M4 ~ 3nT^ one gets from Eq. (42) 

3n{r)T{r)T'{r) = ~iy{r)M3{r) (55) 

as a constitutive equation for the overdamped case. Let us consider separately cases 
with u = 0, and u ^ 0. 

The case: u = 

Then Eq. (|l]) reduces to 
If 1^3 > there is a solution 



Mar = 1^3 — const (56) 



T.(?2^(i + 1,)^" (57) 



V 21 r ro' 

with sonic integration constant vq. The condition 1^3 > means a positive outward 
flow of energy. One thus gets for small r 

2/7 

(58) 



Taking this into account one easily gets from the Poisson equation and hydrostatic 
equilibrium 

n(r) ^ nor-l<'/^ = ^^U, y = 2.2857143, (59) 

1 8 

(/.'(r) ^ cl^or-'/', = -To. (60) 

Notice that density is singular at the origin, but integrable. Thus, this solution can 
be intended to describe a stationary regime without a central mass (as follows from 
Eq. (|60|)), but with a stationary tightening of binaries at the very central region. 
In our setup this last process is reflected through a constant outward heat current. 

The mass insider a sphere of radius r scales as M(r) ^ r^/^. The inverse is 
r ~ M^/s. These are so-called Lagrange radii. If one chooses a set of equidistant 
M values, e.g. M = i/20, i = l,---,20, then the radii will be closer to each 
other near the center than in the outer parts of the cluster. This happen one since 
16/7 > 2. If, on the other hand, the density would be finite at the center, then 
one would have r ~ M^/^, implying large separations between Lagrange radii near 
M = 0. 

The case: u ^ 

Here one can still use Eq. ( |55| ) for M3, which now does not reduce to const/r^. 
Having substituted this equation to Eq. (|4l|), one obtains 

/ = : (r^un)' = 0, 1^1: (nT)' + n(p' + u'un = 0, (61) 

1^2: [r^T^/^T'Y + unT' + 2u'nT = 0, (62) 

7 

For small distances the self-consistent solution of these equations reads 

n(r) ~ r(r) ~ u{r)^uor''/\ (/.'(r) ~ ^Tor-^/^(63) 
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where 6Tq = 57Uono. Notice that there is no central mass (as follows from 
Eq. (p3|)). Moreover, uq has to be positive and u(0) = 0. Therefore, the obtained 
solution can be suitable to describe evaporation phenomena. 



3.3 Underdamped situation 



There is yet another extremal situation realized at sufficiently large distances from 
the center, where the behavior of the system is almost coUisionless (underdamping) . 

The case with u = 

In the underdamping regime one expects that all moments with fc > 3 will be 
equally important, therefore they have to behave in the nearly similar way. Looking 
for a self-consistent, long-distance solution of the hierarchy ( ^9[]4^ ) one gets: 

-3/2 



n(r) 



- j.7/2 ' 



T(r) 



To 



i/(r) ~ 7 



Ma 



V4, v^noT^ 
^ 



-3/2 



(64) 



(65) 



where z/4, no, To are constant. At this stage of our asymptotic analysis they will 
remain unfixed. For 0' one has from the Poisson equation: 

GM SnGno 



4>'{r) 



r5/2 



(66) 



where the first term is the standard asymptotic limit for a finite-mass cluster with 



M as the total mass. The second term follows from Eq. (64). The behavior of M5 
is qualitatively the same as of M4 as seen from Eq. (|42i. Our expression for the 
density given by Eq. (|6j) coincides with that obtained inQ for the halo of a globular 
cluster by means of qualitative physical considerations. 

Having substituted Eqs. (|6^, ^4|) to hydrostatic equilibrium equation (]4C|), one 
gets the correction terms to the density and temperature 



n(r) 



no 
r7/2 



ri/2 



T(r) 



To 
r 



j.l/2 



(67) 



5ro(ni -(- Ti) = GM{ni + SnGno) (68) 

with constant ni, Ti. 

The case with u ^ 

This case is capable to describe the evaporation phenomenon u a as we indicated 
already in section |3.2| . Having this in mind, we will look for a solution, which has 
the Maxwellian structure for very long distances, since there it is supposed to de- 
scribe free escaped stars. In other words, for those scales only first two moments 
u = uo > and T = To > are nonzero, and should be viewed as boundary condi- 
tions. Possible solutions of this kind should be then matched with the overdamped 
solutions having finite u. Notice that the equation (|4^) for energy can be written 
in a simplified form 

^3M' + ^[«'Mnor = o, (69) 



extended'version: submitted to World Scientific on February 1, 2008 11 



where M3 = 1^3 (r) r^, and where J = r^nu is the constant current of particles. 
Searching again for a self-consistent solution, one gets 

2 nio+niilnr 7130 + n-21 Inr + n22(lnr)2 
r n(r) ~ no H \ , (70) 

^/ X ^ , Tio + Tiilnr Tao + Tsi Inr + TaaGnr)^ 

, , 010 + (/'II In r 020 + 021 In r + 022 (In r) 2 
r0~0oH \ ^ , (72) 

Mio+Miilnr M20 + "21 Inr + U22(lnr)2 

u(r)~'UoH \ 5 , (73) 

r 

, . , K310 + J^3ii ln7' 7^320 + J^32i In 1^322 (In r)2 
7/3(r) ~ 1/30 + + , (74) 

.(r)^7^^^^- (75) 

Higher-order terms can be presented analogously. Notice that density does not 
have to be integrable at infinity when considering stationary regimes of evaporation. 
Indeed, the only way to have a truly stationary evaporation regime is to have an 
infinite amount of stars, the most of them being located at infinity. The appearance 
of logarithms in the asymptotic expansions is mathematically connected with non- 
integrable character of the density. 

It is seen as well that the leading-order behavior of vir) is the same for both 
underdamped regimes (with and without evaporation). 

At the present stage of our asymptotic analysis we able to obtain only certain 
relations between the coefficients of ( [70[ - [75| ). 

{r^4>)' = 4:'KGr'^n{r) -> 0o = AnGno, rin = un = 0, 0ii = 47rG'r7,io; (76) 
unr'^ = J, —f J = uqUq, rioUiQ + 7^lo^io = 0, rioMii -|- tihUo = 0; (77) 
(nr)' + 0' + ^(u2)' = O, ^0o = 2ro, (78) 
^0010 - 2710U0U10 = 3rio7io - Sriino -I- 27ZioTo, STn = 0ii; (79) 



J 



u 



[u^T]' = 0, 



J(rii - Tio) - 2T0U10 -I- 7^311 + i^3io = 0, -JTii = 7/311- (80) 

The equations support an additional simplification Tio — 0io = 0, and consequently 
47rGno = Tq — 2uq. This equation connects the energetic characteristics of the 
considered solution. 



4 Conclusion 

This contribution was devoted to the kinetic theory of many-body self-gravitating 
systems in the framework of the model kinetic equation developed by Bhatnager, 
Gross and Krook The main advantage of using this kinetic model (as well as 
other, more refined kinetic models) is that it produces rather tractable schemes in 
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contrast to the full kinetic equations, which are typically quite difficult to inves- 
tigate. On the other hand, model kinetic equations offer rather sensible results, 
which are usually at least in a qualitative agreement with those obtained from the 
full Boltzmann kinetic equation. The physical reason of this success is that a good 
amount of statistical phenomena are quite insensitive to the details of the colli- 
sion integral. In order to explain them it is enough to use the simplest relaxation 
mechanism which only takes into account necessary conservation laws of energy, 
momentum and probability (number of particles) and the correct characteristic 
relaxation time. 

The assumptions on the spherically symmetric character of the considered sys- 
tem and the assumed radial character of non-equilibrium led us to the kinetic equa- 
tion (1), which was then studied by the method of moments. The straightforward 
analytic investigation allowed us to identify several self-similar regimes of a steady 
but non-equilibrium behavior. These solutions corresponds to the presence of the 
central black hole in the system (section 3.1), the influence of binary- formation 
and tightening processes to the system (section 3.2) and evapijration phenomena 
(section 3.3). In this way we reproduced certain known results u and also obtained 
a number of new regimes of behavior. 

We believe that the model kinetic approach will be potentially quite powerful 
in application to many-body, self-gravitating systems. We plan in future to turn 
to the dynamical aspects of this approach, as well as construct more realistic (from 
the viewpoint of applications in astrophysics) kinetic models, that e.g account for 
the strong anisotropy of the velocity distribution in the outer part of the cluster. 
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The kinetic theory of a self-gravitating system is considered in the framework of 
Bhatnager-Gross-Krook model. This approach offers a unique and tractable setup 
for studying the central, collision-dominated region of the system, as well as its 
outer almost coUisionless part. Stationary non-equilibrium states of the system are 
considered and several different self-similar solutions are identified. 



Introduction. There is a loiigptiaxii.tion of applications the kinetic theory to the 
gas of self-gravitating particleafl^oBcl^El. The first regular approach of that kind was 
developed by Chandrasekhar El, and appeared to be very useful in astrophysical 
applications. Indeed, the kinetic theory is certainly applicable for globular clusters 
of stars, which are rather well-isolated systems on the relevant time-scales, and 
contain sufficiently many (typically of order N ~ 10^) stars u. In this contribution 
we present an approach, where one models the coUisional part of the kinetic equa- 
tion, trying to preserve its main qualitative characteristics, and to make the system 
analytically tractable. Our approach is based on one of the most successful models 
,f the kinetic equation proposed by Bhatnager, Gross and Krook (BGK) in fifties 
Since that time it appeared to be very useful, and the reason of this success 
is that the model offers a minimal way to incorporate the necessary conservation 
laws (those of particle number, momentum and energy) with a simple relaxation 
mechanism. This approach provides natural methods to consider simultaneously 
the collision-dominated (overdamped) regions of the system, and its outer almost 
coUisionless domains (underdamping) . In the present contribution we will outline 
the approach, and present some of its results. More details will be found in U. 

Bhatnager-Gross-Krook kinetic equation. The statistical dynamics of N point 
particles with unit mass interacting through Newton's inverse-square law, will be 
considered in the classical kinetic approach. There are two important assumptions: 
(i) The field acting on a given particle (test particle) can be represented as a mean, 
self-consistent field plus a contribution from two-particle collisions, (ii) As usual, a 
collision between a pair of particles is taken to be independent on the others, and 
local (namely, it leads to sudden changes in the corresponding momenta, whereas 
the coordinates can be considered as fixed). In the spirit of the BGK-model we 
propose the following kinetic equation for a spherical self-gravitating system: 



— + -(j) — H 

Ot OVr r 



Vrf 

T{r) 



dl 

dVr 



iy{r) 



n{r)e 



[27rT(r)]i/2 



(1) 



where /(r, Vr) stands for the probability distribution of the radial coordinate r and 
the radial velocity Vr — ^r, is the self-consistent gravitational potential, and 
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where /' = df{r, Vr)/dr, cf)' = d(/)/dr, n{r) — J dvr ,f(r, Vr) and 



n{r)u{r) ^ J dvrVr f{r,Vr), n{r)T{r) ^ J dvr {vr - u) f{r,Vr), (2) 

are the local average radial velocity and the local temperature. Eq. (|l]) is inherently 
related to the Poisson equation: r^^(r^0')' — 47rGn(r). The r.h.s. of Eq. (|l]) is the 
model collision integral, which conserves probability, momentum and energy, and 
enforces relaxation to the corresponding local Maxwellian. ^{r) = ^nT~^^^ (with 
7 cx G^lnA^ as thfipdamping constant) is connected with the inverse characteristic 
relaxation time Elfi^El. Eq. (|l|) is still a complicated integro-diffcrential equation. To 
subtract the relevant information from it one considers the subtracted moments 
Mi{r) = J dwj. [vr — uY f{r,Vr). The corresponding equations read: 



1 / 2r\ 

Ml + liiMi^i + —[r\Mi+i + uMi)]' + I [ (f)' Mi^i + lu'[Mi + uM,_i] 

J.Z \ f J 

= v{r)[u:{l)nT'l\l~iy\-Mi\, (3) 

where iju{1) is zero for I odd, and equal to 1 for I even. 

We will consider the stationary (steady) situation. This case, where all quanti- 
ties are time-independent, is worth to study for its own sake, since there are clear 
observational evidences that typical globular clusters do have such a phase in their 
evolution □'El. On the other hand, it is a necessary step towards understanding of 
more general, time-dependent situation. If the known conditions for equilibrium are 
satisfied (e.g., absence of energy and/or probability currents), the self-gravitating 
system under study has the Gibbs distribution as its unique steady state. As a 
result of the gravothermal instability this distrituition looses its stability for tem- 
peratures lower than certain critical temperature Bfl. However, there are situations, 
where one expects that the very reasons for the existence of the gibbsian equilibrium 
can be invalid. The most typical situation of that kind appears with binary star 
formation and tightening processes in the central region of a star cluster Bu. These 
processes are accompanied with a release of energy, so that the rest of the system 
can be considered as being subjected to sources of energy put in the central part. 
The existence of energy currents from the central part leads to a non-equilibrium 
steady state, and can stabilize the behavior of the self-gravitating system prevent- 
ing its collapse, since the mechanism of the gravothermal instability is connected 
with spontaneous energy currents towards the center 

Ideal hydrodynamics. This scheme is determined by a condition Mj, = (no heat 
current) , and the nonequilibrium character is displayed only through the stationary 
current of particles u. The following exact relations can be inferred from Eqs. (^ 
with Z = 0,1,2: n = cu~^r~'^, T — ciu~^ , where c and ci are some constant 
numbers. If we assume that u does not depend on r, the Poisson equation offers 
an exact solution: (p = 2Tlnr -I- (jiQ. This is an exact non-gibbsian (u ^ 0) but 
isothermal solution having a non-integrable density at the origin. 

The solution with u ^ const has the following form at the small distances: 

u^Uor'l\ <P'^'l^,r-\ T^^r~\ n^^r~^l\ (4) 
2uo Uo "0 
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It is seen that the density is singular at the origin, but stiU integrable. The form of 0' 
indicates the presence of the central mass 7ci/(2MgG'). This implies the following 
choice: mq < 0, c < 0, and then the obtained solution describes a stationary 
consumption regime of stars by the central mass. The energy current nu^ + 3nuT 
is also directed towards the center. Altogether, this solution describes the situation 
with a central mass ("black hole"), but without other dissipative processes. 

Overdamped situation. Starting since this moment, we take into account the 
heat current, described by A/3. We make a self-consistent assumption that there is 
a range of distances close to the inner part (i.e., the core), where collisions dominate. 
In this overdamped regime 77 = (dynamic time) / (kinetic time) is large, and the large 
i'{r) in the r.h.s. of Eqs. (|^) imposes the Gaussian behavior for the corresponding 
moments: M3 = 0(1/ ?y), Ah-SnT^ = 0{l/r]). Using the last resuh as M4 ~ SnT^ 
one gets from Eqs. (^) with I = 3: 

3n{r)T{r)T'{r) = ~iy{r)M3{r) (5) 

as a constitutive equation for the overdamped case. For u = Eq. (|^) reduces to 
Mar^ = 1^3 = const, and for 1/3 > {a. positive outward flow of energy), there is a 
solution which for small r reads 

T^Tor-'/', n(r) ~ nor-i6A, </,'(r) ~ ^or""/', (6) 

where no oc To/(47rG), (/)o esc Tq oc (71^3)^^^. Notice that density is singular at the 
origin, but integrable. Thus, this solution can be intended to describe a stationary 
regime with tightening of binaries at the very central region. In our setup this last 
process is reflected through a constant outward heat current. 

For w 7^ can still use Eq. (|^) for A/3, which now does not reduce to const/r^. 
For small distances the self-consistent solution reads 

n(.)-^, Tir)^^, u{r)^uor-'\ - (7) 

where QT^^^ = 57Mono. Notice that there is no central mass (as follows from Eqs. (||, 
^). Moreover, uq has to be positive and u(0) = 0. Therefore, the obtained solution 
can be suitable to describe evaporation phenomena. 

Underdamped situation. There is yet another extremal situation realized at 
sufficiently long distances from the center, where the behavior of the system is 
close to be collisionless (underdamping). Here one expects that all moments with 
A; > 3 will be equally important, and they have to behave in the nearly similar way. 
Looking for a self-consistent, long-distance solution of the hierarchy (^ with u — 
(no evaporation) one gets: 

t \ "0 rr^f ^ To jno V3 Vi 

- - V' - ' ' " ^' ^' 

where 1^4, no, To are constant. At this stage of our asymptotic analysis they will 
remain unfixed. For (f>' one has from the Poisson equation: 4>'{r) — GMr^^ — 
87rGnor~^/^, where the first term is the standard asymptotic limit for a finite-mass 
cluster with M as the total mass. The behavior of A/5 is qualitatively the same as of 
4. Our expression for the density given by Eq. ^ coincides with that obtained in 
for the halo of a globular cluster by means of qualitative physical considerations. 
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_ Since the case with m 7^ is capable to describe the evaporation phenomenon 
we will look for a solution which has the Maxwellian structure for very large 
distances, because there it describes free escaped stars. For those scales only first 
two moments u = uq > and T = Tq > are nonzero, and should be viewed as 
boundary conditions. Searching for a self-consistent solution, one gets for large r 

rn(r)~no + — — — , r(r) ~ Tq H , 

47rGr r 

r,/)'~2ro + ^i , u{r)^ua^^ , (9) 

r r 

The product r^M^ is a constant to leading order, indicating true energy loss, while 
i>{r) behaves as in Eq. (|8|); Tn and uio are certain constants Notice that density 
does not have to be integrable at infinity, since the only way to get a stationary 
evaporation regime is to have an infinite amount of stars, most of which are located 
at infinity. 

Conclusion. This contribution was devoted to the kinetic theory of many-body 
self-gravitating jSystems in the framework of the Bhatnager-Gross-Krook model ki- 
netic equation In contrast to the full kinetic equation, it produces tractable 
analytical schemes. On the other hand, kinetic models offer sensible results, which 
are usually at least in a qualitative agreement with those obtained from the full 
Boltzmann kinetic equation, since a good amount of statistical phenomena is quite 
insensitive to the details of the collision integral. We report several different self- 
similar regimes of behavior. Some of them are ppw, and in certain cases we confirm 
results obtained by more heuristic approaches □. More details of this investigation 
will be presented elsewhere la. 
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